function [V,T,Mesh] = init_mesh_Xenophontos(p,epsilon,No,Nb,Nc)
%  p----- the method order
% epsilon----then physical parameter
% No ---the number of center element in each row
% Nb ---the number of boundary element in each row
% Nc ---the number of corner element in each row
% gen 1d-mesh distribution , add boundary refinement:
% tao = abs(p*epsilon*log(epsilon));
tao = min(0.25,abs((p+1)*epsilon*log(No)));
% tao = min(0.25,abs((p-1)*epsilon*log(epsilon)));
x1=0:tao/Nb:tao;
x2=tao:(1-2*tao)/No:(1-tao);
x3=(1-tao):tao/Nb:1;
x=[x1(1:end-1) x2 x3(2:end)];

% then generate initial quad mesh
N = length(x);
V = zeros(N*N,2);
% the mesh clumn is : [V1 V2 V3 V4 N1 N2 N3 N4 Active Parent Level]
Mesh = zeros((N-1)^2,11);
% add the init node:
for i=1:N
    for j=1:N
        V((i-1)*N+j,:) = [x(i) x(j)];
    end
end
% and generate the init Quadature mesh
for i=1:N-1
    for j=1:N-1
        temp = zeros(1,11);
        % set the verteics
        temp(1:4) = [(i-1)*N+j (i-1)*N+j+1 i*N+j+1 i*N+j];
        me = (i-1)*(N-1)+j;
        temp(5:8) = [me-N+1 me+1 me+N-1 me-1];
        nei_flag = [1 1 1 1 i-1 N-j-1 N-i-1 j-1];
        pos = find(nei_flag==0);
        temp(pos) = zeros(1,length(pos)); %boundary neighbour is 0
        temp(9) = 1; % set active;
        temp(10) = 0; % the initial elements have no parent;
        temp(11) = 1; % the first level;
        Mesh(me,:) = temp;
    end
end
% here i begin to refine the mesh at the corder
for i = 1:Nc
    [V,Mesh] = refine_corner(V,Mesh);
end
% at last Triangulate it 
Mesh = Quad2Tri(Mesh);
T = get_tri(Mesh);

function [V,Mesh] = refine_corner(V,Mesh)
%find the corner element: only 2 of it's neighbour is 0
rows = find(Mesh(:,9)==1);
flag = or(Mesh(rows,5:8),0);
flag_sum = sum(flag');
corner_element = rows(find(flag_sum==2))'; 
% the row index of corner elements stored in pos
for me = corner_element
    nV = size(V,1); % the new node ranged (nV+,nV+5);    
    %first generate five points:
    V_new = [(V(Mesh(me,1),:)+V(Mesh(me,2),:))/2 ; (V(Mesh(me,2),:)+V(Mesh(me,3),:))/2 ; ...
        (V(Mesh(me,3),:)+V(Mesh(me,4),:))/2 ; (V(Mesh(me,4),:)+V(Mesh(me,1),:))/2 ; ...
        (V(Mesh(me,1),:)+V(Mesh(me,2),:)+V(Mesh(me,3),:)+V(Mesh(me,4),:))/4];
    % and add them to global V
    V = [V ; V_new]; 
    %secondly  add 4 new element:
    nElem = size(Mesh,1);
    % Element data index:  V1  V2  V3  V4   N1  N2  N3  N4  Active Parent Level   
    Mesh = [Mesh; zeros(10,11)]; % one refine step will add 10 elements
    % known info of new elem 1
    Mesh(nElem+1,1:4) = [Mesh(me,1),nV+1,nV+5,nV+4]; 
    Mesh(nElem+1,6:7) = [nElem+2,nElem+4];
    % known info of new elem 2 
    Mesh(nElem+2,1:4) = [nV+1,Mesh(me,2),nV+2,nV+5];
    Mesh(nElem+2,7:8) = [nElem+3,nElem+1];
    % known info of new elem 3
    Mesh(nElem+3,1:4) = [nV+5,nV+2,Mesh(me,3),nV+3];
    Mesh(nElem+3,5) = nElem+2; Mesh(nElem+3,8) = nElem+4;
    % known info of new elem 4
    Mesh(nElem+4,1:4) = [nV+4,nV+5,nV+3,Mesh(me,4)];
    Mesh(nElem+4,5:6) = [nElem+1,nElem+3];
    
    Mesh(nElem+1:nElem+4,9) = ones(4,1);
    Mesh(nElem+1:nElem+4,10) = me*ones(4,1);
    Mesh(nElem+1:nElem+4,11) = (Mesh(me,11)+1)*ones(4,1);
    Mesh(me,9) = 0; %kill this element after new element;
    
    % and treate the neighbour ,there are 2 non zero neighbours
    nei_index = 0;
    %check i's neighbour:
    for i =1:4
      nei = Mesh(me,i+4);
      if(nei~=0)%my first neighbour is not zero then do
        Mesh(nei,9)=0; %kill this neighbour
        elem_begin = nElem+4+nei_index*3;
        Mesh(nElem+i,i+4) = elem_begin+1;
        Mesh(nElem+mod(i,4)+1,i+4) = elem_begin+3; % fill previous unkonws
        nei_index_me = find(Mesh(nei,5:8)==me); % My position in nei_1's neighbour index
        n1 = Mesh(nei,mod(nei_index_me,4)+1);
        n2 = Mesh(nei,mod(nei_index_me+1,4)+1);
        n3 = Mesh(nei,mod(nei_index_me+2,4)+1);
        n4 = Mesh(nei,mod(nei_index_me+3,4)+1);
        %add three trianges and we don't care their neighbour now
        Mesh(elem_begin+1,:) = [n1 n2 nV+i 0 0 0 0 0 1 nei Mesh(nei,11)+1];
        Mesh(elem_begin+2,:) = [n2 n3 nV+i 0 0 0 0 0 1 nei Mesh(nei,11)+1];
        Mesh(elem_begin+3,:) = [n3 n4 nV+i 0 0 0 0 0 1 nei Mesh(nei,11)+1];
        nei_index = nei_index + 1;
      end
    end
end


function Mesh = Quad2Tri(Mesh)
[m,n]=size(Mesh);
for me = 1:m
   if(Mesh(me,9)==1 & Mesh(me,4)~=0)
      % active quadrature
      Mesh(me,9)=0;
      if(length(find(Mesh(me,5:8)==0))==2) %level one
          corner = max(find(Mesh(me,5:8)==0));  %the corner index
          Mesh = [Mesh;Mesh(me,corner) Mesh(me,mod(corner,4)+1) Mesh(me,mod(corner+1,4)+1) 0 0 0 0 0 1 me Mesh(me,11)+1; ...
                   Mesh(me,mod(corner-3,4)+1) Mesh(me,mod(corner-2,4)+1) Mesh(me,corner) 0 0 0 0 0 1 me Mesh(me,11)+1];
      else
          Mesh = [Mesh;Mesh(me,1) Mesh(me,2) Mesh(me,3) 0 0 0 0 0 1 me Mesh(me,11)+1; ...
                   Mesh(me,1) Mesh(me,3) Mesh(me,4) 0 0 0 0 0 1 me Mesh(me,11)+1];
      end
   end
end

function T = get_tri(Mesh)
[m,n]=size(Mesh);
n_tri = length(find((Mesh(:,4)==0)&(Mesh(:,9)==1)));
T = zeros(n_tri,3);
count = 1;
for me = 1:m
   if(Mesh(me,9)==1 & Mesh(me,4)==0)
       T(count,:) = Mesh(me,1:3);
       count = count + 1;
   end
end